# ============================== SETUP ===============================
rm(list = ls())
set.seed(1)
options(scipen = 999)
setwd("~/Dropbox/Wayne-Ying/White_Nationalist_Recruitment/replication/codes")
library(tidyverse)
library(data.table)
library(dplyr)
library(vars)

# ========================= DATA WRANGLING ==========================
comb <- read.csv("../datasets/input/panelbottom25leader.csv", stringsAsFactors = FALSE)

# ============================ VAR MODELS ===========================
# ------------------------ Overall results (follower, leader) -------------------
variables  <- c("follower", "leader", "theme")
mformula   <- as.formula(paste0("~", paste0(variables, collapse = " + ")))
model_data <- model.matrix(mformula, comb[, variables])

# remove intercept
model_data <- model_data[, -1, drop = FALSE]

X_endogenous <- as.matrix(model_data[, c("follower", "leader"), drop = FALSE])
X_exogenous  <- as.matrix(model_data[, setdiff(colnames(model_data), c("follower", "leader")), drop = FALSE])

# Hard filter rows with any NA in endogenous or exogenous
keep <- complete.cases(X_endogenous) & (ncol(X_exogenous) == 0 | complete.cases(X_exogenous))
X_endogenous <- X_endogenous[keep, , drop = FALSE]
X_exogenous  <- if (ncol(X_exogenous) == 0) NULL else X_exogenous[keep, , drop = FALSE]

p <- 2
var_model_merged    <- VAR(y = X_endogenous, p = p, exogen = X_exogenous)
var_irfs_cum_merged <- irf(var_model_merged, n.ahead = 20, cumulative = TRUE)

#===============================================================================
#======================== PART B: EXTRACT THE RESULTS ==========================
#===============================================================================
# Helper to tidy IRF output into wide PE/LWR/UPR
tidy_irf <- function(irf_obj, covar_labels) {
  variables <- names(irf_obj$irf)
  elements_to_pull <- c("irf", "Upper", "Lower")
  out <- NULL
  for (el in elements_to_pull) {
    new_irf_info <- irf_obj[[el]]
    for (outcome in variables) {
      df <- as.data.frame(new_irf_info[[outcome]])
      df_long <- df %>% gather(cov, value)
      df_long$out <- outcome
      df_long$day <- rep(1:nrow(df), length(unique(df_long$cov)))
      df_long$e_type <- el
      out <- rbind(out, df_long)
    }
  }
  out$e_type <- recode(out$e_type, `irf` = "pe", `Lower` = "lwr", `Upper` = "upr")
  colnames(out)[colnames(out) == "value"] <- "val"
  out <- out %>%
    filter(cov != out) %>%
    mutate(val = (val / 10) * 100) %>%
    spread(e_type, val)
  # Relabel for plotting
  out$cov <- recode(out$cov, !!!covar_labels)
  out$out <- recode(out$out, !!!covar_labels)
  out$cov <- factor(out$cov, levels = c("Followers", "Leaders"))
  out$out <- factor(out$out, levels = c("Followers", "Leaders"))
  out
}

# Pooled results (Followers vs Leaders)
irf_data_wide <- tidy_irf(
  var_irfs_cum_merged,
  covar_labels = c("follower" = "Followers", "leader" = "Leaders")
)

#===============================================================================
#========================== PART C: PLOTTING ===================================
#===============================================================================
pdf("../plots/FigureA9.pdf", width = 12, height = 4)

plot_db_pooled <- irf_data_wide %>% dplyr::filter(day == 15)

ggplot(plot_db_pooled, aes(x = cov, y = pe)) +
  geom_segment(aes(x = cov, xend = cov, y = lwr, yend = upr), size = 1) +
  geom_point(size = 3, shape = 16) +
  geom_hline(yintercept = 0, color = "red4") +
  facet_wrap(~ out, nrow = 1) +
  coord_flip() +
  xlab("") +
  scale_y_continuous("\n15-day responses (in percentage points)", expand = c(0.05, 0.1)) +
  theme(
    panel.spacing   = unit(1.05, "lines"),
    legend.position = "none",
    panel.background = element_blank(),
    panel.grid.major = element_line(colour = "gray90", linetype = "solid"),
    axis.text       = element_text(size = 18),
    axis.text.y     = element_text(hjust = 0),
    strip.text      = element_text(size = 20),
    panel.border    = element_rect(colour = "black", fill = FALSE),
    strip.background = element_rect(colour = "black"),
    axis.title      = element_text(size = 16)
  )
dev.off()